{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Parallel stochastic simulation\n",
    "\n",
    "Stochastic simulation is an [embarrassingly parallel problem](https://en.wikipedia.org/wiki/Embarrassingly_parallel) in which realizations are only a function of the random seed, and can be generated indenpendently one from another. Although this is a well-known fact, there has not been reasonable effort in geostatistical software to exploit modern hardware such as HPC clusters and the cloud (e.g. AWS, MS Azure).\n",
    "\n",
    "In GeoStats.jl, *all* stochastic simulation algorithms generate realizations in parallel by default. The package exploits Julia's built-in support for parallel execution, and works seamlessly on personal laptops with multiple cores as well as on high-performance computer clusters with multiple nodes.\n",
    "\n",
    "In this tutorial, we demonstrate how to generate realizations with sequential Gaussian simulation in parallel. The same script can be run on a computer cluster where thousands of processes are available.\n",
    "\n",
    "Before we proceed, please install the following packages:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mPackage GeoStats is already installed\n",
      "\u001b[39m\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mPackage Plots is already installed\n",
      "\u001b[39m\u001b[1m\u001b[36mINFO: \u001b[39m\u001b[22m\u001b[36mPackage PyPlot is already installed\n",
      "\u001b[39m"
     ]
    }
   ],
   "source": [
    "for pkg in [\"GeoStats\", \"Plots\", \"PyPlot\"]\n",
    "    Pkg.add(pkg)\n",
    "end\n",
    "\n",
    "# make sure this tutorial is reproducible\n",
    "srand(2000);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Number of processes\n",
    "\n",
    "When you start Julia, it starts with a single process:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nprocs()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to run simulations in parallel with GeoStats.jl, the first thing we need to do is increase the number of processes in the pool *before* loading the package. The command `addprocs` adds a given number of processes for parallel execution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4-element Array{Int64,1}:\n",
       " 2\n",
       " 3\n",
       " 4\n",
       " 5"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "addprocs()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that `addprocs()` when called without an argument, adds the number of *logical* cores available in the machine.\n",
    "\n",
    "> **WARNING:** If you are on Windows 8 or an older version of the operating system, you will likely experience a slow down. Please add the number of *physical* cores instead that can be found externally in computer settings or via the [Hwloc.jl](https://github.com/JuliaParallel/Hwloc.jl) package.\n",
    "\n",
    "On a HPC cluster, computing resources are generally requested via a resource manager (e.g. SLURM, PBS). In this case, the package [ClusterManagers.jl](https://github.com/JuliaParallel/ClusterManagers.jl) provides variants of the built-in `addprocs()` for adding processes to the pool effortlessly. For example, we can use `addprocs_slurm(1000)` to request 1000 processes in a SLURM job.\n",
    "\n",
    "Now that the processes are available, we can run normal GeoStats.jl scripts and they will automatically distribute the execution among the processes.\n",
    "\n",
    "## Problem definition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2D SimulationProblem (conditional)\n",
       "  data:      3×4 GeoDataFrame (x and y)\n",
       "  domain:    100×100 RegularGrid{Float64,2}\n",
       "  variables: precipitation (Float64)\n",
       "  N° reals:  3"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using GeoStats\n",
    "\n",
    "geodata = readtable(\"data/precipitation.csv\", coordnames=[:x,:y])\n",
    "domain = RegularGrid{Float64}(100,100)\n",
    "problem = SimulationProblem(geodata, domain, :precipitation, 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solving the problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "SeqGaussSim solver\n",
       "  - precipitation => Ordinary Kriging\n"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "solver = SeqGaussSim(\n",
    "    :precipitation => @NT(variogram=SphericalVariogram(range=20.))\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[33mWARNING: \u001b[39m\u001b[22m\u001b[33mSeqGaussSim is not fully implemented, assuming data is already ~ Normal(0,1)\u001b[39m\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "2D SimulationSolution\n",
       "  domain: 100×100 RegularGrid{Float64,2}\n",
       "  variables: precipitation"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "solution = solve(problem, solver)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using Plots; pyplot()\n",
    "\n",
    "plot(solution, size=(1000,300))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.1-pre",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
